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Intrinsic  optical  bistability  in  collections  of  spatially  distributed  two-level  atoms 

Y.  Ben-Aryeh,*  C.  M.  Bowden,  and  J.  C.  Englund 
Research  Directorate,  Research,  Development,  and  Engineering  Center,  Redstone  Arsenal,  Alabama  S5898-S248 

(Received  19  May  1986) 

We  use  a  quantum-electrodynamical,  many-body  treatment  to  show  mirrorless  optical  bistability 
in  terms  of  the  spatial  properties  of  coherent  dipole-dipole  interactions  among  interacting  two-level 
atoms.  The  general  theory  is  applied  to  two  special  cases:  (Da  thin  sample  of  two-level  atoms,  with 
a  width  smaller  than  a  resonance  wavelength  and  (2)  a  long  sample  of  two-level  atoms,  with  dimen¬ 
sions  very  large  relative  to  a  resonance  wavelength.  While  for  the  thin  sample  we  are  able  to  use  a 
mean-field  approximation  with  validity,  for  the  long  sample  we  are  compelled  to  take  into  account 
retardation  and  propagation.  In  both  cases  bistability  is  found  to  be  related  to  a  renormalization  of 
the  frequency  (or  relaxation  rate)  that  is  inversion  dependent.  For  the  long  sample  the  frequency  re¬ 
normalization  is  significant  for  high  atomic  densities  and  for  large  oscillator  strengths. 


-Vi 
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I.  INTRODUCTION 


Intrinsic  optical  bistability  (IOR)  that  is  not  caused  by 
external  feedback  such  as  mirrors,  has  been  the  subject  of 
recent  intense  interest.1  It  was  first  pointed  out  by 
Bowden  and  Sung2"1  and  by  Bowden2"”  that  optical  bista¬ 
bility  (OB)  may  occur  for  a  system  comprised  of  a  collec¬ 
tion  of  atoms  interacting  with  the  electromagnetic  field 
and  driven  by  an  externally  applied  coherent  field  without 
external  feedback.  The  first  detailed  experimental  study 
of  intrinsic  optical  bistability  was  conducted  almost 
simultaneously  by  Hajto  and  Janossy3  using  amorphous 
GeSe2,  who  interpreted  their  results  is  due  to 
temperature-dependent-induced  optical  absorption  in  the 
material,  and  Bohnert,  Kalt,  and  Klingshim,4  who  used 
CdS,  and  also  Rossmann,  Henneberger,  and  Voigt5  using 
the  same  material.  The  process  in  the  first  case  depends 
upon  absorption  due  to  temperature  variation  induced  by 
the  incident  field,  whereas  the  latter  cases  depend  upon  sa¬ 
turation  of  absorption  due  to  the  generation  of  carriers  in 
the  material,  and  the  IOB  has  been  interpreted6  as  due  to 
band-gap  renormalization.  Since  the  earlier  works,  there 
have  been  many  theoretical  and  experimental  investiga¬ 
tions  of  various  forms  of  IOB. 1 

We  present  here  a  fully  quantum-mechanical  treatment 
of  mirrorless  (intrinsic)  optical  bistability  (IOB)  from  a 
collection  of  a  large  number  of  spatially  distributed  two- 
level  atoms  interacting  via  the  electromagnetic  field  and 
driven  by  an  externally  applied  coherent  field.  Atomic 
cooperative  effects  in  IOB  have  been  treated  in  previous 
works  either  by  assuming  a  small  volume  with  dimensions 
smaller  than  a  wavelength,2,7  or  by  assuming  a  system 
with  a  small  number  of  atoms  in  a  semiclassical  approxi¬ 
mation.8  In  the  present  work  we  treat  the  problem  from 
the  many-body  standpoint  by  developing  the  Heisenberg 
equations  of  motion  in  the  “bad-cavity”  limit,  in  which 
the  variables  associated  with  the  field  modes  are  adiabati- 
cally  eliminated. a 

The  general  theory  is  developed  in  Sec.  II.  We  distin¬ 
guish  between  terms  in  the  equations  cf  “s,~>nianeou:,,” 
“cooperative,"  “induced,”  and  “Langevin  force”  origin. 


While  the  spontaneous,  induced,  and  Langevin  force 
terms  depends  only  on  the  single  atom  response,  the 
cooperative  terms  include  factors  of  operators  belonging 
to  pairs  of  atoms  with  arguments  expressive  of  retarda¬ 
tion.  In  the  cooperative  terms  the  explicit  spatial  depen¬ 
dence  of  the  dipole-dipole  interactions  enters  in  the  coeffi¬ 
cients  B(i,j)  which  turn  out  to  be  identical  to  those  ob¬ 
tained  from  classical  dipole-dipole  interactions.10,11 

We  treat  the  steady-state  conditions  for  our  system  and 
by  taking  the  expectation  values  in  the  Heisenberg  equa¬ 
tions  of  motion,  the  Langevin  force  terms  vanish.  In  the 
present  work  we  do  not  treat  fluctuations  and  the  effects 
of  quantum  fluctuations  in  our  system  will  be  presented  in 
a  separate  work.12 

It  has  been  established  in  various  works8,13-15  that  fac¬ 
torization  of  the  cooperative  terms  in  steady-state  results 
in  a  cubic  nonlinearity  and  bistability.  After  establishing 
that  factorization  of  different  atomic  operators  is  a  suit¬ 
able  approximation  for  the  distributed  many-atom  system 
we  show  bistability  effects  under  steady-state  conditions 
for  two  cases:  ( 1 )  the  thin  sample  geometry  with  propaga¬ 
tion  length  smaller  than  a  resonance  wavelength  and  (2) 
the  extended  sample  geometry  where  propagation  effects 
are  important. 

In  Sec.  Ill  we  treat  the  thin  sample  in  the  mean-field 
approximation  in  which  we  use  average  values  for  the  ex¬ 
pectation  values  of  the  atomic  operators.  The  present 
mean-field  approximation  ignores  propagation  effects  and 
this  is  justified  for  a  thin  sample  with  a  propagation  dis¬ 
tance  smaller  than  a  wavelength.  For  the  steady-state 
conditions  a  cubic  equation  is  derived  which  shows  bista¬ 
bility  in  the  parameter  space.  The  linear  stability  analysis 
for  this  system  is  described  in  Sec.  IV. 

In  Sec.  V  we  treat  a  long  sample  by  taking  into  account 
retardation  and  propagation.  We  analyze  the  cooperative 
effects  and  find  that  under  appropriate  approximations 
the  conventional  Maxwell-Bloch  equations  are  reproduced 
but  with  a  correction  that  can  be  expressed  as  an 
inversion-dependent  renormalization  of  the  frequency. 
This  correction  is  proportional  to  the  number  of  two-level 
atoms  per  cubic  wavelength  and  to  the  decay  constant, 
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and  therefore  should  be  important  for  high  densities  and 
strong  oscillator  strengths.  By  taking  into  account  propa¬ 
gation  and  retardation  effects  in  the  long  sample  we  find 
the  interesting  result  that  the  two  phases,  of  high  and  low 
transmissivity,  may  coexist  spatially  in  the  material. 

II.  GENERAL  THEORY 

Our  system  is  composed  of  a  large  number  of  spatially 
distributed  two-level  atoms  coupled  to  each  other  only  via 
the  electromagnetic  field  and  driven  externally  by  an  ap¬ 
plied  radiation  field  taken  to  be  in  a  coherent  state  and 
propagating  along  the  z  axis  with  a  linear  polarization  in 
the  x  direction.  The  Hamiltonian  which  describes  the 
system  is  given  in  the  rotating-wave  and  in  the  electric  di¬ 
pole  approximations  by1 6,17 


H  =H0+H'  , 

N 


Hi 


k  i  =  l 


.V 


—  ift/2  2  ou+e 

i  =  i 


(i)  -<(<“ o'-kori> 


+  H.c.  , 


where  H.c.  denotes  the  Hermitian  conjugate  and 


gl"= 


2iro>k 

fiV 


1/2 

Pp,  x  . 


H0  includes  the  free  atoms  and  free-field  Hamiltonian. 
We  use  here  the  usual  operators  related  to  the  SU(2)  alge¬ 
bra  for  a  single  atom:  Uj"  represents  the  population 
inversion  operator  and  C+1  are  the  raising  and  lowering 
operators  for  atom  /  with  coordinates  r,.  coR  is  the  Rabi 
rate  associated  with  the  applied  coherent  field.  w0  and  k0 
are  the  carrier  frequency  and  wave  vector  of  the  applied 
field,  respectively.  V  is  the  quantization  volume  for  the 
field,  cok  is  the  frequency  of  the  mode  k,  %  is  a  unit  vector 
in  the  direction  of  polarization,  P  is  the  dipole  moment 
matrix  element,  and  p,  is  the  unit  vector  for  the  dipole  of 
atom  i. 

The  Heisenberg  equations  of  motion  for  our  system  are 
obtained  in  the  bad-cavity  limit  by  adiabatically  eliminat¬ 
ing  the  variables  associated  with  the  field  modes.9  We 
eliminate  the  rapid  time  and  spatial  dependence  of  the  di¬ 
pole  operators  by  substituting 


a,;’0(r)  =  ol;)(r)e 


-i<"o'-ko'r,1 


We  get 


(2) 


(3) 


do(J\t) 

dt 


2±' B*(i,j)ouUt)o% 
j  =  < 


f  —  — 
c 


4- o)rCt U')f  +  *(r ,,/)  -f-  H.c. 


d<T+o(t).  =  {ico-02)oi»o(t)+  2'  Bd,j)a% 


dt 


j- 1 


t  — 


<7j\t)e  ~iko'Ui~,s>-L^-sTU> 


+  ~ct,2‘l+/(->(  ij ,  t  )a{‘\  t ) 


The  prime  on  2  indicates  yV=/,  and  /<  +  l(r,,t)  and 
/<_)(r(,t)  are  Langevin  force  operators  which  stem  from 
the  vacuum  contribution  to  the  fluctuations  due  to  normal 
ordering  of  the  field  operators  relative  to  the  atomic 
operators  in  products  of  field  and  atomic  operators  in  the 
equations  of  motion.  These  operators  are  8  correlated, 

{f  +  \rht),?-\rht'))«b(t-l')  .  (5) 

In  our  fully  quantum-mechanical  model  0i=2/9,  P2—P, 
where  f}=2  \  P  \  2k}/3fi  is  defined  as  half  the  spontaneous 
decay  constant.  If  we  introduce  additional  homogeneous 
broadening,  and  (32  may  be  considered  as  empirical  con¬ 
stants.12  In  such  an  approach  the  model  becomes  less 
rigorous,  but  more  general. 

The  coefficients  B(i,j),  derived  from  the  model,  are 
identical  to  those  obtained  from  classical  dipole-dipole  in¬ 
teractions:11 

BUJ) =T0([prPj-(pr?,y  ><p j  lijnwrij ) 

+<P,-VtfrWi<*M  -  (6) 

where 


(4) 

I 

Fi(krij)=e‘kr,)(  1  /k  2rlj  +  i  /k  }r2  -  i  / kr ,  )  , 

(7) 

Fu ( kri}  )  =  elk'ui  —  2 i /k  V,y  —2/k  2r2j )  . 

Here  k  =co/c,  p,  and  p j  are  the  unit  vectors  for  the  di¬ 
poles  of  the  atoms  i  and  j,  rl;  =  r,  —  Tj  and  r,;  =  |  r,  —  r;  | 
is  the  distance  between  the  two  atoms. 

In  developing  Eqs.  (3)  and  (4)  we  sum  the  interaction  of 
atom  i  with  many  other  atoms  described  by  the  summa¬ 
tion  over  j  ( j=£i).  In  the  present  work  we  assume  a  model 
based  upon  many-body  interactions  in  which  each  atom  is 
influenced  by  the  ensemble  average  over  all  other  atoms. 
Such  an  assumption  is  justified  when  the  number  of  atoms 
within  a  volume  of  a  wavelength  is  large  (  » 1 ).  We  ap¬ 
ply  the  equations  either  for  steady-state  conditions  or  for 
small  fluctuations  near  the  steady  state  (i.e.,  linear  stabili¬ 
ty  analysis).  Under  such  conditions  any  correlation  which 
is  generated  between  the  atoms  due  to  initial  conditions 
(like  that  assumed  in  the  transient  behavior  of  superradi- 
ance)  is  destroyed  by  the  dephasing  mechanism.  The  time 
of  decorrelation  is  of  order  1  /&  where  P  is  the  spontane¬ 
ous  decay  time  and  becomes  shorter  if  we  add  additional 
homogeneous  broadening  due  to  collisions.  Under  such 
conditions,  as  shown  recently  by  Hopf  and  Bowden8  in 
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numerical  simulations  for  a  finite  number  of  atoms  within  <  a>|»0( ,  )a</>0(  t-rij/c))  =  (  a%(  r)><  a%(  i  -r,  /c))  , 

a  volume  of  a  wavelength,  we  can  adequately  justify  fac-  >J  ’  (8) 

torization  of  the  products  of  the  dipole  operators  among  /  /  ui,  ,  lW  w,  , , 

different  atoms.  Factorization  is  further  Sidated  by  the  -nj/c)az  U»  =  W0U -r./clXa^U))  , 

fact  that  there  is  a  distribution  of  retardation  times.  etc.  By  taking  nto  account  that  expectation  values  of 

'W  .  _  1*  .  .  . .  .1  .  _  .•  v  ■  »  .  _  _  .  _ 


Therefore  we  use  the  approximation 


Langevin  force  terms  vanish,  [Eq.  (4)],  we  get 


d(az(t)){ 


=  -A[<tri,,(r)>  +  l]-  <uR(a%(t))  +22'  B'UJ)el^ir,-'>\a%(t)){o%U  -r0/c))  +c.c. 


c/( -%(()) 


=  (iA— /?2)<cr'{o(t)>  +(<uj /2)<ffj  (r)>  +  ^'B{i,j)e 


;^-,ko'(r(-ry' 


(a't'oU  ~riJ/c))(a'z‘,U))  , 


where  c.c.  denotes  the  complex  conjugate  and  A=<u— <y0 
is  the  deviation  of  the  applied  field  frequency  from  reso¬ 
nance.  The  reaction  field  is  defined  as 

©(r,)=  -r^/c))  .  (11) 

7  =  1 

The  cooperative  terms  in  Eqs.  (9)  and  (10)  represent  the 
effect  of  the  reaction  field,  which  for  low  values  of  the 
externally  applied  field,  reduces  the  internal  field.  For 
higher  values  of  the  external  field  its  value  decreases  sud¬ 
denly  and  thus,  as  shown  below,  produces  a  first-order 
phase  transition.  In  order  to  illuminate  the  physical  prop¬ 
erties  of  our  system  we  apply  our  general  equations  to  two 
special  cases:  (Da  thin  sample  of  two-level  atoms,  with  a 
width  smaller  than  a  resonance  wavelength  X  and  (2)  a 
long  sample  of  two-level  atoms,  with  dimensions  very 
large  relative  to  a  wavelength. 

III.  THIN  SAMPLE 

For  coherent  radiation  impinging  on  a  thin  film  of 
two-level  atoms  with  a  large  surface  area  and  a  width 
smaller  than  a  wavelength  X  we  ignore  the  time  of  retar¬ 
dation  rtj/c  by  a  Taylor  expansion  of  the  operators  and 
use  the  approximation 

.  (12) 

Under  these  approximations  we  can  use  mean  values  for 
the  expectation  values  of  the  atomic  operators: 

<<7+'o (t  — '■y/c)>~(<7%(f)>~<cr%(7)>  =  <«7  +  0(/)>  , 

< <r'i o( t  — fjy/c))~<olio(f) )~<<7(lo(t))  =  <er_0(r)>  . 

After  using  Eqs.  (12)  and  (13)  in  Eqs.  (9)  and  (10),  we  get 

d(crz ) 

— —  =  -0,(<az>  +  l)-4Rer<<7+o><ff_o> 

-«R<(7+o>-Wr<CT_o)  •  (14) 

— -^^-=(iA—f}1)(cr+0)  +  r(a+0)(oz)  +  ~-(oz )  , 


where  all  expectation  values  are  taken  at  the  same  time  t 


r=  2'  b  uj). 

7  =  1 

The  mean-field  approximation  [Eq.  (13)]  becomes  a  very 
good  approximation  when  T  has  only  a  small  dependence 
on  the  location  of  each  atom  i  in  the  thin  sample.  In  such 
cases  one  can  use  Eqs.  (14)  and  (15)  with  some  mean  value 
for  T.  In  order  to  check  the  validity  of  this  approxima¬ 
tion  and  also  in  order  to  give  some  estimates  for  the  value 
of  T  we  have  made  the  following  calculation. 

We  assume  that  the  thin  sample  is  a  cylindrical  slab 
with  a  thickness  d  «A..  The  z  direction  is  defined  as  the 
direction  k  of  propagation  of  the  external  field  and  the  di¬ 
poles  induced  by  the  externally  applied  field  are  assumed 
to  be  oriented  along  the  jL  axis.  As  described  in  Appendix 
A  we  calculate  T  as  a  function  of  z,  where  (0,0,  z, )  is  the 
location  of  the  atom  i  and  where  the  origin  of  the  coordi¬ 
nates  is  assumed  to  be  at  the  center  of  the  cylinder.  The 
general  result  represented  in  (A6)  is  quite  complicated  but 
for  the  limit  kd-*0,  ke-+0,  d/e=const,  we  get  a  simple 
analytical  result: 


—  3ifinX} 


(d2/4-z,2)1/2 


—  d/2+e<Zj  <d /2  —  e  , 

where  n  is  the  density  of  the  two-level  atoms.  In  the  cal¬ 
culation  we  excluded  a  volume  y7T63=l/n  about  r,  from 
the  integration  range  in  order  to  exclude  the  self-field  of 
atom  i.  T  depends  on  z,  only  through  the  logarithmic 
function  which  is  a  slowly  varying  function.  This  con¬ 
clusion  is  quite  good  even  for  atoms  which  are  very  near 
the  surface  [based  on  the  comparison  between  Eqs.  (A7) 
and  (A8)]. 

For  the  thin  film  the  main  contribution  to  the  reaction 
field  comes  from  the  dipole  field  in  the  “near  region.”10 
If  the  system  has  dimensions  larger  than  a  wavelength 
and  if  the  dipoie-dipole  interaction  is  of  spherical  or  cubic 
symmetry  the  average  contribution  of  the  dipole-dipole  in¬ 
teraction  to  the  near  region  vanishes.  This  is  similar  to 
the  condition  in  the  usual  derivation  of  the  Lorentz- 
Lorenz  correction.10  However,  for  a  thin  film  with  a 
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width  much  smaller  than  a  wavelength  this  local  symme¬ 
try  is  broken  due  to  the  width  of  the  film,  and  the  contri¬ 
bution  of  the  near  region  field  becomes  dominant.  We 
find  that  in  this  case  the  value  of  T  diverges  unless  we  ex¬ 
clude  the  self-field  from  the  volume  of  integration. 

Under  the  steady-state  conditions  we  get  from  Eqs.  (14) 
and  (15): 


<o+0>  =  77 


-(.(ol/  2){<jz) 


+u  («A-ft)+r<af>  ’ 

0,(  < a, )  +  1 )[( < az  >Rer-02)2  +  ( < az ) Imr  +  A )2] 

=  —02  I  I  2<ffi  >  ■  (18) 

Equation  (18)  is  a  cubic  equation  in  the  population  inver¬ 
sion  and  as  proved  in  the  next  section,  it  leads  to  bistabili¬ 
ty  effects.  The  critical  points  for  (a2)  as  a  function  of 
|  co/j  | 2  are  determined  by  taking  the  derivative  of  Eq.  ( 1 8) 
with  respect  to  (crz)  and  setting  d  |  o)R  1 2/d(az  >  =0. 
We  get 

<ffz>2+F(<7z>+G=0,  (19) 

-thrift  U  (20, 

3  in2  3  in2 

t(A2+^)-2Re(r)ft  +  2Im(r)A  +  (ft//3,)|wR  |2] 


This  must  be  satisfied  simultaneously  with  Eq.  (18).  The 
condition  for  the  threshold  is  G  =  jF2  and  by  assuming 
the  approximations  r  »/32,  T  »  |  A  | ,  we  get 

Alfiil.i  (22) 

@i  I  r | !  J 

for  the  threshold  condition. 

We  find  that,  as  the  field  intensity  increases  and  the 
Rabi  parameter  becomes  large  relative  to  the  cooperative 
constant  T,  the  two-level  system  switches  suddenly  from 
the  low-transmission  branch  to  the  high-transmission 
branch.  We  find  here  the  interesting  result  that  the 
switching  occurs  at  relatively  lower  intensities  if  dephas¬ 
ing  mechanisms  due  to  collisions  are  introduced.  This  is 
illustrated  in  the  next  section  in  Fig.  4. 

In  previous  works, l,2,7,s  it  was  shown  that  nonlinear  re¬ 
normalization  of  the  frequency  may  lead  to  bistability.  In 
the  present  case  such  renormalization  is  obtained  from 


Eqs.  (14)  and  (15)  in  which  A— »A  +  (az  )Im(D. 

According  to  the  theoretical  calculations  made  in  Ap¬ 
pendix  A,  T  becomes  imaginary  in  the  limit  that  kd  — 0, 
d/e= const.  By  performing  numerical  calculations  of  T 
for  a  thin  sample,  we  found  that  the  real  part  of  T  was 
small  relative  to  the  imaginary  part  and  that  the  ratio  de¬ 
creases  as  a  function  of  the  density  of  atoms.  While  the 
imaginary  part  of  T  leads  to  a  renormalization  of  the  fre¬ 
quency,  the  real  part  of  T,  which  is  positive,  is  related  to 
superradiance.  According  to  the  present  analysis  both  the 
real  part  and  imaginary  part  of  T  may  lead  to  bistability. 
In  treating  the  bistability  effect  for  the  thin  sample  one 
should  note  that  (o+0)  is  proportional  to  cj*r  and  it  will 
tend  to  zero  only  for  a  vanishing  external  field  amplitude. 

IV.  LINEAR  STABILITY  ANALYSIS 
FOR  THE  THIN  SAMPLE 

Here  we  show  that  the  solutions  to  the  steady-state  con¬ 
ditions  of  Eq.  (18)  indeed  exhibit  bistability.  Let  us  con¬ 
sider  a  stationary  state  of  the  thin  sample.  For  infini¬ 
tesimal  perturbations  of  the  system  from  the  stationary 
state  the  linearized  set  of  equations  obtained  from  Eqs. 
(14)  and  (15)  are 


db(az ) 


=  -/3,6(  crz  }  -4Ren6(cr+o>)<o_o> 

—  4Rer(o+0>(6(o_o>  > 

—  wR(b(cr  +q))~ a)Rb(a _q)  , 


db(a+0) 


=  (/A— &)8<  a+o  >  +  r(8<<7+0>  )(crz  ) 

4-  r(cr+0}(5(<72 ))  +  («« /2)5(<7Z )  ,  (24) 


where  (o+0)  and  (az)  are  again  given  by  Eqs.  (17)  and 
(18).  The  equation  for  5(o_0)  is  the  complex  conjugate 
of  Eq.  (24). 

By  considering  the  three-component  vector  q 
=  (6(crz),8(a+0),6(o'_o))  and  following  a  general  pro¬ 
cedure  for  linear  stability  analysis,18  we  introduce  in  the 
linearized  equations  the  ansatz 

q( /)  =  q(0)exp( Ar)  ,  (25) 

where  A  is  a  complex  number.  We  get  the  polynomial 
equation  for  the  stability  eigenvalues, 

A.3-MiA.2-M2A.-M3=0,  (26) 


i —  Pi  +  2/32  —  2 Re( r)(<7j )  , 

^2  =  (0i/&)(  <<** )  +  1  )Im(D[A  +  Im(  n((72 )  ]  +  (3/?i//?2)Re(  D(  (oz }  +  1  )[Re(  D((7Z )  —02] 

+  [(<trz)Rer  — /32)2  +  ((oz  )ImT  +  A)2][  1  —  (/3t//?2)((<7z  >  +  1  )/(<tz  )] 

+  2^1ft-2^Re(r)<aI>-2(ReD2(l3,/ft)<CTl)(<az>-|-l)  , 

^3  =  2#l((oj )  +  1 )[  |  T  |  2(ct2  )  — /32Re(r)  + A  Im(r )]—(/?! /(ct2  > ) j ( (trz  jRef  — )32)2  + [ (oz  )Im( D-t- A]2 1 
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In  deriving  Eqs.  (27)  we  eliminated  (<x+0),  (a_o)>  and 
|  a)R  | 2  by  the  use  of  Eqs.  (17)  and  (18). 

According  to  the  Hurwitz  criterion19  the  polynomial  of 
Eq.  (26)  has  only  zeros  with  negative  real  parts  if  the  three 
conditions 

^,>0,  A,A2-A 3>0,  /»3>0  (28) 

are  fulfilled.  These  are  therefore  the  conditions  for  our 
system  to  be  stable  under  small  shifts  from  a  stationary 
state. 

Since,  according  to  Eq.  (18),  the  equation  for  (az)  at 
the  turning  points  is 


-A3 

Mo,)'0’ 


(29) 


the  condition  of  stability  .43  > 0  implies  that  the  system 
is,  as  expected,18  unstable  in  the  region  in  which  we  have  a 
negative  slope: 


d(az)/d  \oR  |2<0  . 


For  positive  values  of  Ret  I")  the  condition  A\  >0  is  al¬ 
ways  fulfilled  (as  <cr2 )  is,  in  our  system,  always  negative). 
Our  system  becomes  unstable  in  regions  of  positive  slope20 
when  Re(D  is  negative  and  2Re(D(<r2)  Al¬ 

though  in  our  theoretical  and  numerical  calculations, 
Re(D  was  positive,  there  might  be  cases  of  subradiance 
for  which  Re(  T )  is  negative  and  the  instabilities  obtained 
in  such  cases  should  be  of  interest. 

In  order  to  illustrate  the  properties  of  the  steady  states 
in  the  thin  sample  we  have  calculated  (az)  as  a  function 
of  | 2  [Eq.  (18)]  and  the  stability  conditions  [Eqs. 
(28)]  for  various  examples  in  which  we  vary  the  physical 
parameters.  Some  representative  examples  are  pictured  in 
Figs.  1—5  where  the  solid  and  dotted  curves  represent, 
respectively,  stable  and  unstable  steady  states.  We  have 
taken  0,  =2,  so  that  all  rates  are  essentially  in  units  of 
(/?t/2);  thus  |  coR  | 2  is  in  units  of  (/?,/ 2)2. 

While  in  Figs.  1—4,  Re(D>0  and  the  instabilities  are 


(°i 


FIG.  2.  ( az )  as  a  function  of  |  a>R  | 2  for  a  thin  sample  with 
normalized  parameters:  4?,  =2,  Pi—\,  Re(D=0,  Im(D=  — 10. 
The  deviation  from  resonance  is  given  by  (a)  A  =  —  8,  (b) 
A=-6,  (c)  A=-4,  (d)  A=-2,  (e)  A=-l,  (f)  A=0,  (g)  A  =  j, 
(h)  A=l.  The  solid  and  dotted  curves  represent,  respectively, 
stable  and  unstable  states. 


only  in  the  regions  of  a  negative  slope,  in  Fig.  5,  Re(D  <0 
and  the  instabilities  are  also  in  regions  of  a  positive 
slope.20  These  figures  demonstrate  therefore  the  condi¬ 
tions  A  |  >  0  and  A  3  >  0  for  the  stability  of  our  system. 
The  condition  A  j  A  2  >  A  3  was  fulfilled  in  all  the  numeri¬ 
cal  calculations  in  which  ,4 ,  >  0  and  A  3  >  0  and  therefore 
did  not  give  additional  instabilities. 

According  to  the  theoretical  analysis  the  figures  are  in¬ 
variant  to  simultaneous  changes  in  the  signs  of  Im(  T )  and 
A.  As  demonstrated  in  these  figures  the  bistability  is  re¬ 
lated  to  the  imaginary  and/or  the  real  part  of  T. 

In  Figs.  1—3  and  5  we  assumed  the  relation  /?2=/?i/2 
which  follows  from  the  fully  quantum-mechanical  model 
in  which  the  spontaneous  decay  stems  from  the  self-field 


0  M  |  i-  log  1M 

FIG.  1.  (<r2)  as  a  function  of  |  £l>r  | 2  for  a  thin  sample  with  1,1 


normalized  parameters:  0,=2,  ft=I,  A=0,  Re(D=0.  The  FIG.  3.  (<72 )  as  a  function  of  | o)r  | 2  for  a  thin  sample  wiih 
imaginary  value  of  T  is  given  by  (a)  T/= —4,  (b)  T/ = —6,  (c)  normalized  parameters:  /5]  =  2,  02=1,  A=0,  lm(D=  —  10. 
T/  =  —  8,  (d)  T;  =  -  10,  (e)  T,  =  —  12,  (f)  T/  =  —  14.  The  solid  The  real  part  of  T  is  given  by  (a)  T*  =0,  (b)  T*  =2,  (c)  T*  =4, 
and  dotted  curves  represent,  respectively,  stable  and  unstable  (d)  T*  =6,  (e)  T*=8,  (f)  T*  =  10.  The  solid  and  dotted  curves 
states.  represent,  respectively,  stable  and  unstable  states. 
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FIG.  4.  ( az )  as  a  function  of  |  ojr  | 2  for  a  thin  sample  with 
normalized  parameters:  f}\=2,  A  =0,  Re(D=0,  Im(D=  — 10. 
The  dephasing  constant  is  given  by  (a)  /3,  =  2,  (b)  /32=  1 .5,  (c) 

A=l. 


part  of  the  atom-field  interaction.  For  taking  into  ac¬ 
count  additional  homogeneous  broadening,  one  may  con¬ 
sider  the  ratio  /32//3|  as  an  empirical  parameter  and  the  ef¬ 
fect  of  changing  this  parameter  is  illustrated  in  Fig.  4. 

The  threshold  condition  given  in  Eq.  (22)  is  in  fairly 
good  agreement  with  our  numerical  calculations.  The  ex¬ 
plicit  threshold  condition  for  bistability  in  terms  of  Im(  F ) 
shown  in  Fig.  1  is  consistent  with  similar  threshold  condi¬ 
tions  determined  in  Refs.  7  and  8.  Further  conclusions  on 
the  qualitative  behavior  of  our  system  can  be  obtained  by 
examining  the  systematic  dependence  of  the  figures  on  the 
various  system  parameters. 


V.  LONG  SAMPLE  WITH  RETARDATION 
AND  PROPAGATION  EFFECTS 


For  coherent  radiation  transmitted  through  a  long  sam¬ 
ple,  with  dimensions  very  large  relative  to  a  wavelength, 
we  use  again  Eqs.  (9)  and  (10)  and  study  the  effect  of  the 
spatially  varying  parameters  on  the  behavior  of  the  sys¬ 
tem.  We  define  the  z  axis  as  the  direction  of  propagation 
of  the  externally  applied  field,  and  the  x  axis  as  the  direc¬ 
tion  of  polarization.  Since  the  dipoles  induced  by  the 
externally  applied  field  are  approximately  in  the  x  direc¬ 
tion  we  simplify  the  expression  given  in  Eq.  (6)  for  B(i,j) 
by  assuming 


p,  =  p;  =  Px  .  (30) 

For  calculating  the  cooperative  interaction  of  atom  i 
with  all  other  atoms  we  separate  this  interaction  into  two 
parts.  We  choose  a  sphere  around  the  atom  i  with  a  ra¬ 
dius  r0  that  is  on  the  order  of  a  wavelength,  so  that  within 
this  sphere  we  can  ignore  retardation  and  use  the  mean- 
field  approximation.  The  second  part  of  the  cooperative 
_ l 


_  3  0nk3 

y  =  R  -  jsin(2 R)- 


1  +  cos(2R)  2sin(2R) 

R  R2 


1  —  cos(  2  R ) 

R3 


t 


! 


FIG.  5.  (ox )  as  a  function  of  |  o>/5  | 2  for  a  thin  sample  with 
normalized  parameters:  /3,  =  2,  /?2  =  1,  A  =  0,  Im(F)=  — 10. 
The  real  part  of  T  is  given  by  (a)  T*  =  —  2,  (b)  r*  =  —  6,  (c) 
r*  =  -10.  The  solid  and  dotted  curves  represent,  respectively, 
stable  and  unstable  states. 


interaction  includes  the  interaction  with  all  other  atoms 
that  are  not  located  inside  the  sphere.  For  the  second  part 
of  the  interaction  we  take  into  account  retardation  and 
propagation  effects  but  simplify  the  calculations  by  using 
the  approximation  r,;  >  A.,  where  rf;  is  the  distance  be¬ 
tween  atom  i  and  any  other  atom  j  that  is  outside  the 
sphere.  We  transform  the  summations  to  integrations  by 
using  spatially  continuous  variables  for  the  atomic  expec¬ 
tation  values.  We  calculate  the  total  cooperative  effect  by 
adding  the  two  parts  of  the  interaction  and  the  result,  as 
expected,  is  found  to  be  independent  of  our  arbitrary 
choice  of  the  radius  r0  of  the  sphere. 

By  using  the  mean-field  approximation  of  Eq.  (13)  we 
express  the  first  part  of  the  cooperative  interaction  as 

[d(a{J\t))/dt]coopU)  =  -4(a%U))(a%(t))Rtrl  , 

[d { o"U t))/dt ]coop,( ,,  =  < a'j'ot t)){ax\t))V |  , 
where  for  ko~k  =  kz  we  get 

r,=  21'>B(,-,y')e-'*u'-V.  (32) 

j 

In  Appendix  B  we  calculate  T ,  for  the  atom  ;  which  is  at 
the  center  of  a  sphere  with  a  radius  r0  and  exchange  the 
summation  by  integration  for  randomly  distributed  two- 
level  atoms  with  a  density  n.  In  principle,  one  should 
eliminate  the  self-field  by  excluding  from  the  integral  a 
small  sphere  around  atom  i  with  a  radius  e  where 
(47r/3)e  n  =  1.  Since  the  integrand  is  well  behaved  for 
small  r  (see  Appendix  B),  we  assume  the  approximation 
e  =  0  for  ke«  1. 

The  final  result  for  T |  given  by  Eq.  (B4)  can  be  written 

as 


(33) 


—  i 


-7  +  jcos(2R)- 


sin(2R) 

R 


2  cos(  2 R )  sin(2R) 

R2  +  R3 


34 


INTRINSIC  OPTICAL  BISTABILITY  IN  COLLECTIONS  OF  .  .  . 


3923 


where  R  =kr0. 

The  real  and  imaginary  parts  of  y  are  described  in  Fig. 
6  as  functions  of  r0/k.  Ignoring  small  oscillations  we 
find  that  Re(D  is  given  approximately  by  a  linearly  in¬ 
creasing  function  of  r0, 

Re(y)~-^-r0=/{  .  (35) 

As  we  increase  the  radius  r0  of  the  sphere,  the  contribu¬ 
tion  from  the  far  dipoles  becomes  more  dominant  and  the 
effect  of  the  oscillatory  part  of  Re(y )  becomes  negligible. 

Im(y)  is  a  negative  oscillating  function  of  r0/k  with  an 
averaged  value 

Im(y)=  —  j  ,  (36) 

which  is  independent  of  r0  for  r0/X>  1.  The  imaginary 
part  of  y  is  therefore  contributed  by  the  dipoles  in  the  “in¬ 
termediate  region”  and  its  effect  on  the  equations  of 
motion  is  to  produce,  through  Eq.  (31),  a  renormalization 
of  the  resonance  frequency.  According  to  Eq.  (33),  the 
frequency  renormalization  is  proportional  to  the  density 
of  the  two-level  atoms  and  to  the  oscillator  strength. 

For  calculating  the  second  part  of  interaction  with  the 
distant  dipoles  which  are  outside  of  the  sphere,  we  should 
take  into  account  retardation  and  propagation  effects. 
For  this  purpose  it  is  convenient,  at  this  stage  of  the  cal¬ 
culation,  to  describe  the  expectation  values  of  the  atomic 
operators  as  continuous  variables.  We  define  <ff+0(r,t)) 
and  (<7j(r,f)),  respectively,  as  the  complex  dipole  and  the 
inversion  of  population  per  unit  volume  at  point  r,  con¬ 
sistent  with  the  volume  of  integration  of  the  sphere  of  ra¬ 
dius  r0.  Equations  (9)  and  (10)  can  be  written  in  the  new 
variables  as 


FIG.  6.  Re(y)  and  Im(y)  described  as  functions  of  r  /k. 


For  calculating  the  second  part  of  the  reaction  field,  we 
use  the  “radiation-zone  approximation”10  for  B(i,j), 
which  is  valid  for  r y  >  X: 

3  R 

B  ( i,j)  ~  -  - - sin2a  ,  (40) 

2  krjj 


where  a  is  the  angle  between  r,;  and  the  x  axis.  The 
second  part  of  0  can  be  written  as 


e2(r,t)=  — 


3i P 


J'dVj(7  +  o 


r ',t  - 


r-r' 


-'V"— r- >  ,*  ,r  — r*  | 

* - rT~ — m - sin2a  .  (41) 

k  r-r' 


d{az(r,t)) 

It 


+n] 


-  j[20*(r,t)-t  bjR  ]<a+0(r,t))  +c.c.  j  , 


(37) 


d(a+0(r,t)) 

dt 


=  (iA-/32Xcr+o(r,f)> 


+ 


0(r,f)  + 


coR 


(<7j(r,/)>  , 


(38) 


where  0(r,t)  is  a  continuous  variable  corresponding  to  the 
definition  of  reaction  field  0(r,,t)  given  in  Eq.  (11),  and  is 
separated  into  two  parts.  0  =  0, +02  corresponding  to 
the  contributions  within  and  without  the  sphere.  Accord¬ 
ing  to  Eq.  (31),  we  get 


0i(r,t)  =  <o+o(r,r))r|/n 


(39) 


The  prime  on  the  integral  indicates  that  we  have  to  ex¬ 
clude  from  the  volume  of  integration  the  sphere  which  in¬ 
cludes  the  first  part  of  the  interaction. 

For  calculating  02(r),  we  assume  that  the  external ly  ap¬ 
plied  radiation  is  nearly  resonant  with  the  two-level  atom¬ 
ic  frequency  so  that  k  ~&0,  and  find  that,  because  of  the 
spatial  phase  factors  the  main  contribution  to  the  integral 
comes  from  the  region  for  which  k0‘(r  —  r')~  f  r — r'  ( . 
The  strength  of  the  coherent  dipole-dipole  interaction  for 
each  atom  is  therefore  mainly  due  to  atoms  that  are  locat¬ 
ed  at  smaller  values  of  the  propagation  distance  in  the 
material.  Following  these  considerations,  we  assume  that 
the  expectation  values  of  the  atomic  operators  per  unit 
volume  depend  on  the  z  coordinate  and  use  in  Eq.  (41)  the 
approximation  sin2a=  1,  |  r— r'  ]  ~  | z  —z'  \  (for  the 

denominator)  and  |  y  —y'  \  «  |  z  — z'  | ,  |  x  —x'  \ 

«  |  z  —  z'  |  (for  the  exponents).  Following  these  approxi¬ 
mations  we  get,  by  a  straightforward  integration  of  Eq. 
(41 )  over  the  x  and  y  coordinates,  the  result 


02(r,f)  =  02(z,f)  = 


-3Q3  r2- 
2  Jo 

r‘-'o 
k 2  Jo 


dz 


dz’( 


r 


<cr+0(z'+  — (r  -z')/c)> 
k(z-z') 

+0(z',/  —  (z  -z')/c))  . 


f  f  dx  dy  exp\ik[(y  —y')2  +  {x  -x')2]/2(z  -z')\ 


(42) 
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Wc  find  according  to  Eq.  (42)  that  for  each  atom  the 
dipole-dipole  interaction  with  other  atoms  decreases  as  the 
inverse  of  distance  but  that  due  to  phase  factors  appearing 
in  the  equations,  the  number  of  atoms  in  the  opposite 
direction  to  that  of  propagation  of  the  external  field, 
which  contribute  to  the  cooperative  effect,  increases  pro¬ 
portional  to  distance.  So  the  explicit  dependence  of  the 
cooperative  dipole-dipole  interaction  on  distance  is  com¬ 
pletely  eliminated. 

The  choice  of  r0  defines  the  spherical  volume  over 
which  0i  is  calculated.  This  volume  must  be  larger  than 
a  cubic  wavelength  in  order  that  the  radiation-zone  ap¬ 
proximation  be  valid  in  the  calculation  for  02.  On  the 
other  hand,  it  should  be  small  enough  so  that  retardation 
effects  are  not  important  within  this  sphere.  By  using  Eq. 
(36)  we  find  that  the  average  value  of  Im(y)  is  indepen¬ 
dent  of  r  for  r/A.>  1.  The  effect  of  Im(y),  which  enters 
in  Eq.  (3d)  via  Eq.  (33),  is  very  important  as  it  leads  to  a 
renormalization  of  the  frequency.  Following  these  con¬ 
siderations  we  define  r0  to  be  of  order  X,  but  somewhat 
larger  than  X 

Concerning  Re(y )  one  should  note  that  the  effect  of 
Re<  y )  in  0,  is  overwhelmed  by  the  large  value  of  02  in 
the  long  sample,  since  both  terms  influence  the  relaxation. 
Also  by  ignoring  retardation  in  the  small  volume  of  radius 
rQ  and  by  neglecting  the  small  oscillations  of  Ref  y )  in  Eq. 
(34)  we  find,  according  to  Eqs.  (31)  and  (39),  that  the  net 
effect  of  Refy)  is  to  allow  the  upper  limit  of  integration 
for  02  to  become  z  instead  of  z  —  r0. 

We  arrive  now,  after  the  above  calculations,  to  an  im¬ 
portant  physical  result.  Since  (cr+0(z',t  —(z  —z')/c))  of 
Eq.  (42)  is  taken  at  retarded  time  the  derivative  d02/dz 
gives  the  form  of  Maxwell’s  equation  in  the  slowly  vary¬ 
ing  envelope  approximation  and  in  the  retarded  time 
frame.  However,  because  of  the  macroscopic  self-field 
contribution  0|(r)  we  should  take  into  account  in  Eq.  (38) 
frequency  and  dephasing  renormalization  that  is  propor¬ 
tional  to  the  population  inversion  (az( r)).  This  result  is 
remarkable  since  in  all  textbooks  on  quantum  optics  in 
which  Maxwell-Bloch  equations  are  used,  this  effect  is 
neglected.  The  present  analysis  shows  that  it  may  be  im¬ 
portant  for  high  densities  and  large  oscillator  strengths. 

In  summary,  the  equation  for  02(z,/)  is  given  by  the 
form  of  the  ordinary  Maxwell  equation  in  the  slowly 
varying  envelope  approximation: 


I  ae2(*'f)  dQzUf) 

C  df  +  0Z 


~{o  +o(^,f)  )  , 


while  in  the  Bloch  equations  represented  by  Eqs.  (37)  and 
(38)  one  should  substitute  0  =  0|+©2,  where  01  leads  to 
a  renormalization  of  the  frequency,  relaxation,  and  de¬ 
phasing. 

In  the  present  theory  the  spherical  volume  over  which 
the  discrete  variables  are  averaged  to  obtain  the  continu¬ 
ous  variables  can  be  of  order  rw  >  X.  The  use  of  the  fac¬ 
torization  and  mean-field  approximations  for  calculating 
0i  is  justified  only  for  densities  which  are  sufficiently 
large  so  that  we  have  many  atoms  within  a  volume  of  a 
wavelength. 

In  previous  works,1, 2,7,8  it  has  been  shown  that  a  renor¬ 
malization  of  frequency  can  lead  to  bistability.  In  the 


present  work  we  have  clarified  the  nature  of  this  effect  in 
a  many-atom  system  and  found  that  the  renormalization 
of  frequency  is  proportional  to  the  decay  constant  fi  mul¬ 
tiplied  by  the  number  of  two-level  atoms  per  cubic  volume 
of  wavelength,  consistent  with  the  results  of  Refs.  7  and  8. 

If  one  introduces  in  Eqs.  (37)  and  (38)  the  definitions 

(Or  (  T )  =  (Or  -}-  202  (  T)  , 

,  ,  (44) 

r'  =  ©,(r)/<(7+o(r)  > 

~  — i (3/iTr2)-^Xip  , 

these  equations  are  of  the  same  form  as  Eqs.  (14)  and  (15) 
for  the  thin  sample,  with  an  imaginary  value  for  T.  The 
threshold  condition  for  bistability  represented  in  Eq.  (22) 
will  be  correct  here  also  when  we  exchange  o>r  and  T  by 
co'r  and  T',  respectively.  In  the  long  sample,  |o))i(r)|  will 
be  a  decreasing  function  of  the  distance  along  the  z  axis. 
We  get  therefore  the  important  result  that  the  two  phases 
of  high  and  low  transmissivity  may  coexist  spatially  in  the 
material. 

VI.  CONCLUSIONS 

In  the  present  work  we  have  presented  a  general 
quantum-mechanical  treatment  of  coherent  dipole-dipole 
interactions  in  a  coherently  driven  collection  of  a  large 
number  of  spatially  extended,  two-level  atoms,  and  have 
shown  that  a  frequency  renormalization  (as  well  as  relaxa¬ 
tion  renormalization)  that  is  proportional  to  the  popula¬ 
tion  inversion  exists  and  may  lead  to  intrinsic  optical  bi¬ 
stability.  In  the  context  of  the  general  theory  we  clarified 
by  our  analysis  the  relation  between  these  effects  and  the 
use  of  the  conventional  Maxwell-Bloch  equations. 

We  conclude  from  the  present  analysis  that  intrinsic 
(mirrorless)  optical  bistability  may  be  produced  by  a  thin 
sample  of  two-level  atoms  with  a  width  smaller  than  a 
wavelength,  due  to  a  iocal  spatial  symmetry  breaking.  A 
system  of  a  thin  sample  with  two-level  atoms  might  be 
realized  in  experiments  with  Rydberg  atoms.  While  the 
present  analysis  is  suitable  mainly  for  an  atomic  system,  it 
suggests  that  similar  effects  may  be  obtained  by  using 
crystals  with  a  high  density  of  bound  excitons  with  large 
oscillator  strengths.  The  treatment  of  the  long  sample 
shows  the  possibility  that  two  phases,  of  high  and  low 
transmissivity,  may  coexist  spatially  in  the  material. 

According  to  the  analysis  in  the  present  article,  the  ef¬ 
fects  of  intrinsic  bistability  follow  mainly  from  the 
inversion-dependent  frequency  renormalization.  For  the 
long  sample,  this  conclusion  follows  from  Eqs.  (37)  and 
(38),  with  the  value  of  0,  =  r'<o+o)  given  by  Eq.  (44) 
(where  T'  is  imaginary).  For  the  thin  sample,  it  follows 
from  Eqs.  (17)  and  (18),  with  the  value  of  F  estimated  in 
Eq.  (A7)  (the  effect  of  relaxation  renormalization  ex¬ 
pressed  by  the  real  part  of  T  is  negligible  for  high  densi¬ 
ties).  In  order  to  get  me  phase  transition  in  the  two  cases, 
the  frequency  of  the  coherent  radiation  should  be  some¬ 
what  lower  than  the  atomic  frequency  and  of  such  a  value 
that  the  frequency  renormalization  will  bring  it  in  and  out 
of  resonance  as  a  function  of  the  value  of  (az).  To  ob¬ 
tain  a  good  condition  for  intrinsic  bistability,  ImT  should 
be  made  an  order  of  magnitude  (or  more)  larger  than  /?2, 
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where  ft  's  the  relaxation  rate  for  the  complex  dipole. 
Assuming  for  the  pure  quantum  case  ft  =  ft  where  13  is 
one-half  of  the  spontaneous  decay  rate,  we  find  that  we 
need  approximately  50—100  (or  more)  atoms  within  the 
volume  of  a  cubic  wavelength  in  order  that  the  effects  will 
be  significant.  If  we  introduce  additional  homogeneous 
broadening,  which  increases  the  value  of  ft,  the  density 
should  be  increased  by  an  additional  factor  of  ft//?.  The 
present  effects  should  therefore  be  significant  for  high 
densities  and/or  large  oscillator  strengths  and  when  the 
mechanism  of  spontaneous  decay  is  strong.  The  intensity 
of  the  external  field  for  which  the  phase  transition  occurs 
can  be  estimated  by  (22). 

In  the  present  work  we  have  used  the  mean-field  ap¬ 
proximation  in  the  small  volume  to  evaluate  0,,  Eq.  (39), 
which  enables  us  to  treat  intrinsic  bistability  by  an  analyt¬ 
ical  quantum-mechanical  calculation.  The  validity  of  the 
present  model  is  supported  by  previous  treatments  of  this 
problem  from  other  points  of  view.  The  use  of  a  semi- 
classical  model7  consisting  of  the  Maxwell-Bloch  formula¬ 
tion  for  a  medium  consisting  of  laser-driven  two-level 
atoms  with  the  local-field  correction,  where  the  mean- 
field  and  slowly  varying  envelope  approximations  were 
not  made,  gives  precisely  the  same  results  for  the  renor¬ 
malization  of  frequency  and  the  same  threshold  condi¬ 


r=yftiJo  d<t>  J  d/2dz  J ^  pdp[sin2<t>  Fx(kr)  +  cos2<b  Fn(kr)]  ,  (Al) 

where  r  =[p2  +  (z  —  z,)2]l/2.  In  addition,  we  exclude  a  volume  (4/3)7re3=  1/n  about  r,  from  the  integration  range.  Then, 
performing  the  angle  integration,  we  find 


tions  for  bistability  in  the  thin  sample  as  the  present  work. 
The  results  of  the  heuristic  model8  which  treats  the  intrin¬ 
sic  bistability  and  the  dynamical  statistical  behavior  by 
numerical  simulations,  for  a  finite  number  of  atoms 
within  a  cubic  wavelength  quantitatively  confirms  the 
essential  deterministic  results  derived  in  the  present  work, 
as  well  as  that  of  Ref.  7.  However,  the  present  treatment 
is  more  general  than  the  previous  work,  and  is  also  analyt¬ 
ical  and  therefore  should  be  of  much  interest.  We  hope 
that  experiments  will  be  done  for  observing  the  phenome¬ 
na  predicted  theoretically  by  the  present  analysis. 
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APPENDIX  A:  CALCULATION  OF  T  FOR  THE  THIN 
SAMPLE 

The  present  theoretical  calculation  is  done  by  using  the 
approximation  (13)  and  assuming  ko~k  =  A:z.  Locate 
atom  i  at  (0,0, z, )  and  assume  cylindrical  symmetry,  where 
the  origin  of  coordinates  is  in  the  center  of  the  cylinder. 
Then 
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Now  assume  kd-+  0,/ce-*0,d/e  =  const.  Then 

r=  —  2i  -t-/  [i2  — — — .)./  [*'*  — — —  =. 
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for  d/2  +  e<zi  <d/ 2-e.  For  |  z,  |  =  c(/2  we  get  by  a  similar  calculation,  the  result 


r-W"“ 


APPENDIX  B:  CALCULATION  OF  T  FOR  A  SPHERICAL  SAMPLE 


Let  r,  =0  and  use  spherical  coordinates.  Then,  with  ko=)t  =  kz,  and  by  using  the  approximation  of  Eq.  113),  we  have 


rt=  B(i,j)e,kz,=  j0n  fQ  d<t>  Jq  dd  dr  r2sin9e‘krco'e[(s\n2dsm2(!>  +  co‘i2d)F\(kr)  +  sin‘dcos2<l>Fn(kr)] 
=  \n^n  fQ  d6 sin#  J rm‘*  dr  r2sin2de'krc‘yie[(sin2&  +  2cos2d)Fl(kr)  t-sin2VFn(krt]  . 


Let  x  =  cos  <9;  then, 


T,  =  y Trfin  J  "'“  dr  j~  dx  r2cos{krx)\Fi(kr)  +  Fll(kr)  +  x2[Fi(kr)  —  Flllkr)]\ 


=  3 ir/3n  f  dr  r2  (sinkr)/kr[Fl(kr)  +  Fu(kr)]+  —  CPS^-  -  r-~  -sinAir  [F^kr)  — Fu(kr)] 
J  'mm  Jfc  rl  ksri 


or,  with  kr  —R, 


r,=  6 rr£n  j^m“d/{|sin/j[1/*  —  3 //? 3  —  /  ( 1  —  2/7? 2  +  3//?4)]  +  cos.R  [3/R 2  — /  ( 1 /R  -l/Rl)]\e,R 
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Since  the  integrand  is  well  behaved  for  small  R  [it  is  O ( /J 2 ) -4- *0 (/?)],  we  take  Rmin=0  and  let  Rmix—>R.  The  integra¬ 
tion  then  gives 


r,  3 fink2  „  i  .  l+cos(2R)  2sin2R  1  —  cost 2/? ) 

1|  = - 7—  R  -Tsin(2R)  — - - - + - - — - - r - 
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,  .  7  i  sin( 2/2 >  2 cost 2 R)  sin( 2/? ) 
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